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Abstract 


The feasibility of adding viscoelasticity and the Generalized Method of Cells 
(GMC) for micro-mechanical viscoelastic behavior into the commercial 
HyperSizer structural analysis and optimization code was investigated. The 
viscoelasticity methodology was developed in four steps. First, a simplified 
algorithm was devised to test the iterative time stepping method for simple 1-D 
multiple ply structures. Second, GMC code was made into a callable subroutine 
and incorporated into the 1 -D code to test the accuracy and usability of the code. 
Third, the viscoelastic time-stepping and iterative scheme was incorporated into 
HyperSizer for homogeneous, isotropic viscoelastic materials. Finally, the GMC 
was included in a version of HyperSizer. MS Windows executable files 
implementing each of these steps is delivered with this report, as well as source 
code. The findings of this research are that both viscoelasticity and GMC are 
feasible and valuable additions to HyperSizer and that the door is open for more 
advanced non-linear capability, such as viscoplasticity. 
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Summary 

During the course of this work, the feasibility of including viscoelasticity and the 
Generalized Method of Cells (GMC) into a commercial version of HyperSizer was 
demonstrated through two simplified codes and a special research version of HyperSizer. 
There were no prohibitive method, analytical, or software roadblocks to implementing 
viscoelasticity and GMC into the commercial version of HyperSizer. Further, we see a 
potential for future research that could extend the present methodologies to 
viscoplasticity or non-linear post buckling analyses. 

When the process of incorporating the GMC was first begun, the GMC code was setup in 
such a way as to store large matrices within the routine for each analyzed ply. This led to 
very large memory requirements (40 MB per ply). This problem was resolved by the 
subcontractor. Dr Orozco, by solving each of these matrices on the fly and not storing 
them. The final version was able to solve for 100 plies using only 23 megabytes of 
memory. In addition to problems with memory, we were concerned with the amount of 
time that the process would take as a large matrix inversion was being done every time a 
new ply was sent to the GMC routine. This problem was solved using a sub-iteration 
technique that allowed single plies to be analyzed multiple times and prevent the matrix 
inversion from occurring for each call to the GMC routine. 

Demonstrated Technologies 

First, with the development of the Generalized Method of Cells (GMC) code spanning 
most of 1 997, an algorithm was designed to implement this technology into the 
commercially available code called HyperSizer. The algorithm, outlined in the 
attachment, "Progress on Laminate/Panel Viscoelastic Integration," was applied first to a 
simplified problem involving viscoelasticity in 1-D material plies with homogeneous 
isothermal viscoelastic and elastic materials. This simplified algorithm was developed 
into a standalone MS Windows (95/98/NT) executable called Multiplies. This program 
is delivered with this report as a Microsoft Visual Basic source code and Windows 
executable file. Multiplies accomplished several things: 

■ The feasibility of the iterative, time stepping algorithms depicted in the in the 
attachment, "Progress on Laminate/Panel Viscoelastic Integration." 

■ The feasibility of a method devised by Collier Research Corporation and the 
subcontractor. University of Virginia to increase the efficiency of the basic forward- 
Euler time stepping scheme. This method involves doing sub-iterations at each time 
step to increase the accuracy without adversely affecting the total run time. 

■ In the course of development, several closed form benchmark solutions were 
developed against which some of the later HyperSizer generated results could be 
compared. 

Second, the GMC FORTRAN code delivered by the sub-contractor was distilled into a 
callable subroutine for inclusion into the HyperSizer main program. Before incorporating 
this code into HyperSizer, the callable routine, called EpsEtaGMC was incorporated into 
the simplified code, Multiplies to arrive at a new code called MultipliesGMC. This 
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code, also a Visual Basic source code with Windows executable is delivered with this 
final report. MultipliesGMC demonstrated the following: 

■ The feasibility of the iterative, time-stepping algorithm with sub-iteration technique 
applied to the simplified viscoelastic problem using the UVa. generated GMC method 

■ The usability of the GMC callable subroutine for a separate program. 

■ When compared with the closed form solutions and the solutions generated with 
Multiplies, the accuracy of the GMC methodology was demonstrated, at least for 
homogeneous material problems. 

Third , the HyperSizer program was modified to add a viscoelastic capability using the 
same algorithm as in Multiplies, albeit a 3-D version. The first tests were performed 
without GMC to simplify the procedure. The major task involved here was to add a time- 
stepping ability. As HyperSizer normally is restricted to doing point-analyses, it 
previously had no ability to step through time. In the course of this implementation, 
attention must be paid to the histories of the viscoelastic plies. Because HyperSizer 
analysis is based on analysis objects, and each object carries unique strain fields and is 
composed of may plies, die stresses, strains and viscous stresses and strains is tracked for 
every ply as part of every analysis object. To put this in perspective, if a complicated 
panel concept, such as a hat-stiffened panel is composed of composite panels with, say, 

18 plies for the facesheet and six plies for the hats, the number of separate plies tracked 
as part of the analysis is over 60. That is, 60 plies must be calculated for every panel. 

This analysis capability demonstrates that: 

■ The iterative time-stepping algorithm with sub-iteration technique is feasible in a 
HyperSizer full panel implementation. 

■ Again using closed-form benchmarks, the HyperSizer implementation of the Visco- 
elastic solution demonstrates consistent and stable results 

Finally, the HyperSizer program with viscoelastic capability was modified to include the 
same callable subroutine module EpsEtaGMC that implemented the UVA GMC 
FORTRAN code. This modification was relatively minor based on the earlier work 
performed to implement viscoelasticity into the main HyperSizer code and the testing and 
research that went into building and testing the GMC routine through MultipliesGMC. 
This final code, Hs_Visco, is delivered with this final report as a MS Windows 
Executable. The GMC FORTRAN subroutine, EpsEtaGMC and two other supporting 
routines are included as electronic files. This code really demonstrates the promise of 
including the GMC methodology into the HyperSizer commercial code. Although the 
code delivered with the final report is still considered to be 'research' code, it does tell us 
that there are no roadblocks to implementing viscoelasticity and GMC into the 
commercial version of HyperSizer. 


Commercializing HyperSizer Viscoelasticity 

Implementing viscoelasticity in HyperSizer proved to be valuable even without the GMC 
methodology. It is feasible to solve problems in which only certain parts of a complex 
panel would be viscoelastic in nature. For example, for a panel only the facesheet might 
go viscoelastic while the stiffeners remained elastic. The consistent stiffness 
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formulations of HyperSizer ensure that during this process, the forces would remain in 
balance over the entire panel. At the present time, there are two shortcomings in the 
viscoelastic formulation that would have to be addressed for commercialization. 

First, the current implementation of the viscoelastic parameters requires that each ply of 
each analysis object be entered as a separate material. While the histories of each of 
these plies would certainly need to be kept, there is no reason that multiple plies could 
not point to a single material. As implemented right now, the material for each ply is 
entered through a manually generated ASCII data file and the amount of data that must 
be generated for a complex panel concept is prohibitive. For an entire vehicle with many 
individual panels, the amount of data that would have to be generated for many panels 
would be impractical. If multiple plies pointed at a single material, the amount of data 
needed would be more manageable. 

Second, in order to be commerciable, the viscoelastic capability would have to be 
integrated with the current database /GUI scheme. This would require generation of 
viscoelastic material data to be integrated into the existing material editor forms and 
time-step/total time/sub-iteration control to be added to the existing project editor forms. 
In addition to integrating the viscoelastic parameters into the interface, support would 
have to be added to allow for material properties that vary with temperature. Currently 
all viscoelasticity in HyperSizer is isothermal. 


Commercializing HyperSizer Viscoelasticity with GMC 

The Generalized Method of Cells methodology created by the sub-contractor has proven 
to be valuable in its ability to model the micro-mechanical behavior of any general uni- 
directional composite material. The commercialization potential for this capability exists 
but is a little farther away. There are several things that must be accomplished with the 
GMC to make it a part of the commercial HyperSizer code. 

First, just as with the homogeneous viscoelastic algorithm, the material properties for the 
current implementation of the GMC must be entered for every ply of every analysis 
object. This information must be simplified such that multiple plies can point at the same 
material and the materials data must be integrated into the HyperSizer database/graphical 
user interface scheme. This is complicated by the fact that the current GMC subroutine 
reads in its own data file (separate from HyperSizer) and keeps track of all of its own data 
in the course of an analysis. This is inconsistent with the way data structures in 
HyperSizer, which are database driven, are handled. Therefore, the input of all material 
data must be transferred from the GMC routine to the HyperSizer GUI. In addition, the 
GMC routine opens several output data files and writes data into them. Again, this is 
inconsistent with the HyperSizer implementation and the GMC must be modified to give 
all input and output control to the HyperSizer calling progr am 

Second, in addition to just the FilelO, it would be preferable if the control of the history 
data and all allocable arrays were turned over to HyperSizer. Presently, the GMC 
routine, while having a clean interface with HyperSizer, keeps all of its history and ply 
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data separately from HyperSizer. It would be more efficient and cleaner (and easier to 
modify) if the data structures were kept track of by HyperSizer and passed to the GMC 
routine as needed. 

Third, the GMC routine is currently limited to 100 plies. This means that the code would 
only handle one or two panels at the same time and is the biggest challenge that would 
have to be overcome to make the code truly commerciable as part of HyperSizer. One of 
the strengths of HyperSizer is the ability to model a complex structure, such as a 
complete airframe, very quickly. The GMC limitation of 100 plies would be a major 
limitation of this ability and must be overcome for successful commercialization. 


Overall Commercialization Potential 

We see the commercialization potential of the viscoelastic and GMC capability to be 
high. As discussed above, this would require additional effort on our end in the 
integration and interface to make the capability more tightly coupled with the way 
HyperSizer works. In addition, to make the GMC commerciable would require effort on 
the part of the sub-contractor to make the GMC routine fit more tightly into the 
HyperSizer database schema by turning over control of the material and history data to 
HyperSizer and increasing the number of plies. Other than these implementation issues, 
there were no major roadblocks found to adding a commerciable viscoelastic capability 
into HyperSizer. 


Feasibility 

The effort of this contract has definitely demonstrated the potential for using 
viscoelasticity and GMC in HyperSizer. The algorithms developed and implemented for 
viscoelasticity in HyperSizer have demonstrated that the coarse-mesh geometry modeling 
techniques used by HyperSizer are compatible with both viscoelasticity and the GMC 
method. 

This research has also demonstrated that HyperSizer is compatible without limitation 
with non-linear, iterative, and time stepping schemes such as viscoelasticity. Using 
similar techniques there are not roadblocks for extending this capability beyond the 
current viscoelastic methodology into viscoplasticity. In addition to viscoplasticity, 
other non-linear methods could also be implemented using these techniques. For 
example, there has recently been interest by HyperSizer users in post-buckling strength 
analysis for panels and beams. We see the present research as providing a powerful base 
from which to launch such capabilities. 
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Viscoelasticity Implementation in HyperSizer 

Using the some time-stepping algorithm as outlined in the attachment, "Progress on 
Laminate/Panel Viscoelastic Integration," the capability to perform viscoelastic 
calculations on individual plies was added to HyperSizer. HyperSizer uses the concept of 
'analysis objects' to carry out its thermo-mechanical failure and stress analyses. An 
analysis object might be the clear span between hat stiffeners for a hat stiffened panel, the 
web of a hat stiffener, the hat crown, the closed span of the facesheet over the hat, or the 
combo object composed of the flange and facesheet combined. For a composite 
material, each of these analysis objects is made up of multiple plies, and for the 
viscoelastic analysis, the history of each ply in each analysis object must be kept track of. 
For the homogeneous viscoelastic material analysis a data structure and pointer system 
was created in HyperSizer to hold this history and implement a forward Euler time- 
stepping scheme. 

In order to use this methodology, a separate MS Windows executable file is included 
with this report called Hs_Visco.exe. This file is located in the Program 
Files\HyperSizer\Executable directory. (Note: This file is not part of the normal 
HyperSizer installation). The procedure for setting up a viscoelastic run in HyperSizer is 
very similar to setting up a regular HyperSizer run. First, the user sets up a normal 
HyperSizer analysis run (see the HyperSizer documentation for details) and runs the 
normal HyperSizer analysis module to initialize the analysis. Next, a file is created in the 
project materials directory called *. VMT (Viscoelastic MaTerial) which has viscoelastic 
properties for every ply in the analysis. To date, since these files are being completed by 
hand, only simple single panel analyses have been performed. An example * VMT file 
for a three-ply analysis is listed here: 


8000. Total Tima (s) 

200. Tima Step siza (s) 

1 Number of sub-iterations 

0 GMC analysis flag ( 1 - GMC; 0 ■ no 0*2) 

Ply Number (Hot Epoxy) 

1 

Damping Coefficient; ViscoElastic Stiffness; Elastic Stiffness 

58013050., 446701., 48586. 

Ply Number (Cold Epoxy) 

2 

Damping Coefficient; ViscoElastic Stiffness; Elastic Stiffness 

21319797., 588397., 28861. 

Ply Number (Hot Epoxy) 

3 

Damping Coefficient; ViscoElastic Stiffness; Elastic Stiffness 

58013050., 446701., 48586. 

ENDFILE 


The beginning of the file contains the viscoelastic analysis control parameters, total 
analysis time, time step, the number of sub-iterations to perform per time step and a 
control over whether the GMC is to be performed. For each ply to be analyzed, the 
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(lamping coefficient, r| (psi^s 1 ), the elastic stiffness coefficient, E s (psi) and the viscous 
stiffness coefficient, Em (psi). Note that these stiffnesses override the stiffnesses that are 
specified in the HyperSizer Graphical User interface. This is a disconnect in the 
integration of the viscoelasticity into HyperSizer and will be resolved in any follow-on 
development. Once this file is setup, simply go to the Project Input directory and drag 
the Project.ini file onto the executable, Hs_Visco to run the viscoelastic analysis. At the 
end of analysis, an object based plot of stress distribution will automatically pop up to 
indicate the stress distribution among the various plies. 


Results 

This code was used to solve a simple problem involving an unstiffened plate composed of 
three homogeneous viscoelastic plies with uniform compression loading. This is 
implemented as a simple panel concept from the unstiffened plate/sandwich family with 
the * .VMT file from the previous page. In order to compare the results of this analysis to 
those presented in the attached report, "Progress on Laminate/Panel ViscoElastic 
Integration," a transverse load, N y = -v*N x was applied so that the response of the panel 
would be solely in the X (longitudinal) direction. This was done using the user-defined 
load capability in HyperSizer. 



N x = -2400 lb/in 


Running the analysis results in the following plot that shows the amount of force in each 
of the three objects: 


Aqogoo n nnnnn” ,<iril ™ riri nn "™ rifinr *" nwrmnrinnnnnnc 




Obpct 1 
Otwd 2 
Qbm* 3 
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The example performed here has two viscoelastic materials. The center ply, which is less 
viscoelastic (i.e. Em is lower) is twice as thick as the outer plies. This makes sure that 
the stack is symmetric and that the stack and will not experience any bending loads 
during the analysis. As the figure illustrates, during the analysis, creep occurs and the 
less viscoelastic center ply picks up more of the total force while the outer plies pick up 
less force. The summed force (indicated as "Object") remains constant as it should in an 
analysis of creep. The same case run with the standalone code. Multiplies results in the 
following normalized stress distribution: 



Note that this matches almost exactly with the plot given for the HyperSizer analysis. 
Since this is a creep test, it is useful to plot the solution for no rmalize d strain over the 
course of the analysis. 
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These results match very closely with the normalized strains achieved with the 
Multiplies code: 



The results obtained with this non-GMC implementation of viscoelasticity demonstrate 
that the capability has been successfully integrated into HyperSizer, at least for very 
simple geometry problems. The methodology has been demonstrated for more complex 
panel shapes such as the integral blade stiffened panel. For more complex shapes, while 
there are no exact solutions against which to compare, it is worth noting that the initial 
deformed shape of the panel can be determined using the total elasticity of each 
viscoelastic material (Etotai = Em + E s ). The end point (asymptotic) condition of the panel 
can easily be determined by using the standalone elastic portion of the material stiffness, 
Es. Using these material properties to obtain the beginning and ending conditions of the 
panel, and ensuring that at all times the total forces in each panel object always sum to 
the total panel force give reassurance that the method is integrated properly. 

While this method is integrated and gives good time dependent results for viscoelastic 
analysis, it is still a good distance from being a commerciable component of HyperSizer. 
First of all, there is no integration with the HyperSizer database or the graphical user 
interface. This component is essential for a commercial code to let the user avoid having 
to generate ASCII data files by hand in which it is easy to make a typing mistake and 
there is no data integrity assurance. Second, material properties at this time must be 
entered for each ply of each analysis object of the panel concept For the simple case 
illustrated here, there are only three analysis objects consisting of one ply each. For a 
more complicated geometry, many more ply definitions would be required. For example, 
a hat stiffened panel, with an 18 ply facesheet and 6 ply hat stiffener requires over 60 ply 
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material definitions. That is 60 definitions for a single panel! If there are hundreds of 
panels across a vehicle surface, the amount of input definition data becomes prohibitive. 
Therefore, more integration and efficiency work is required to make this implementation 
commerciable. 

GMC Viscoelasticity Implementation in HyperSizer 

With the basic time-stepping and iteration algorithm in place, a callable GMC routine 
called EpsEtaGMC was formed from the GMCVisc program written by Dr.Corlos 
Orozco. This routine was almost a separate program on its own which read its own input 
data from a data file, kept track of its own history data and basically returned the viscous 
strain terms (s' 1 ) given the current strains and previous time step viscous strains. The 
technique proved to work well, though the separation of the routine data from that of the 
main HyperSizer program could prove to make full integration into the commercial 
version more difficult. 


To perform the GMC based viscoelastic analysis, the procedure is the same as the 
previous viscoelastic setup except that an additional data file must be setup for the GMC 
input data. Again, the viscoelastic material properties for each ply was entered into the 
* .VMT file, except that the GMC analysis flag was this time entered as T instead of 'O'. 
In addition, the input file for the GMC analysis looked like this: 

4 0 Mx«Db«r of materials, viacoplst. flag | 

'glass-v 

2 Number of different temperatures 

1 0.2100D+02 Temperature (1) 

0 . 3447D+05 0.3447D+05 0.2872D+20 E«,Em,Eta 

0 . 2000D+00 0 . 2000D+00 Nuaxial , Nutrans 

0.1000D-05 0 . 1000D-05 alphaxial, alp h trans 

2 0 . 1210D+03 Teeperature (1) 

0 . 3447D+05 0.3447D+05 0.2872D+20 Es,Em,£ta 

0 . 2000D+00 0 . 2000D+00 Nuaxial, Nutrans 

0 . 1000D-05 0 . 1000D-05 alphaxial, alphtrans 

r epoxy- r {cold) f 

2 Number of different temperatures 

1 0 . 2100D+02 Temperature (i) 

0 . 40S7D+04 0 . 1990D+03 0.1470D+06 Es,ttn,Eta 

0.3110D+00 0.3110D+00 Nuaxial , Nutrans 

0.1000D-05 0.1000D-05 al phaxial , alphtrans 

2 0 . 12100D+03 Temperature (i) 

0 . 3080D+04 0 . 3350D+03 0.4000D+06 Es,Em,Eta 

0 . 3170D+00 0.3170D+00 Nuaxial , Nutrans 

0 . 10000-05 0 . 1000D-Q5 alphaxial, alphtrans 

* epoxy-r (hot) * 

1 Nxsnber of different temperatures 

1 0 . 2100D+02 Temperature (i) 

0 . 3080D+04 0 . 3350XH03 0.4000D+06 Es,Em,Eta 

0 . 3170D+00 0 . 3170D+00 Nuaxial , Nutrans 

0.1000D~05 0 . 1000D-05 alphaxial, alphtrans 

' epoxy-r (hot) 2 ' 

1 Number of different teaperatures 

1 0.2100D+02 Temperature (i) 

0 . 3080D+04 0 . 3350D+03 1.000D+20 Ss,Em,Eta 

0.3170D+00 0.3170D+00 Nuaxial , Nutrans 

0.1000D-05 0.1000D-05 alphaxial, alphtrans 


Initialization Data and 
material definitions 
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Number of FLY 

0 . 5000D-02 0 . 1000D-05 

2 100 
2 


strain amp, strrate 
100 loadopt , nint , nintact 

50 nstep, natepa 


0 .2100D+02 

4 


0 . 1000D+01 raftanp, daltatamp 

4 ncallay, ncallar 


3 

3 

3 

3 

10 


Geometry ID 


0.0020 
0.0020 
0.0020 
0.0020 
0 .5000D+01 


2 2 2 2 
2 2 2 2 
2 2 2 2 
2 2 2 2 

10 

0.0020 

0.0020 

0.0020 

0.0020 

0.5000D+01 

4 

3 3 3 3 
3 3 3 3 
3 3 3 3 
3 3 3 3 

10 

0.0125 
0.0125 
0.0125 
0.0125 
0 . 5000D+01 


0 . 5000D+01 0 . 5000D+01 0.5000D+01 

4 ncellay , ncellax 


Geometry ID 


0 . 5000D+01 0 . 5000D+01 0.5000D+01 

ncellay, ncellax 


Geometry ID 


0 . 5000D+01 0 . 5000D+01 0.5000D+01 



This file, located in the project Material directory is called <Project>.GMC. Note that the 
number of plies defined in this GMC file must match the number of plies defined in the 
*.VMT file For more information about the GMC file, see the GMC Vise documentation 
listed in Appendix B. 


Note that there is currently a disconnect in the material stiffness calculation between the 
* VMT file and the * .GMC file. This disconnect is in the afact that the material 
stiffnesses are not passed back from the GMC module and therefore the only stiffnesses 
that HyperSizer knows about are those defined in the VMT file. For homogeneous 


N AS A/CR— 200 1-211166 


11 



viscoelastic materials, this is not a problem. However, for composite plies with two or 
more dissimilar materials, there is no way for HyperSizer to know about the stiffnesses 
that the GMC module is using. If this research effort is continued, this deficiency must 
be addressed. 


Results 

In order to test the implementation of the GMC viscoelasticity, a problem identical to the 
three-ply model from the previous section was setup with HyperSizer using the GMC 
methodology. The same *.VMT file from the previous problem was used (with the GMC 
analysis flag set to 1 ) and the *.GMC file as shown on the previous page was used. The 
run parameters and material properties were identical. The stress plots from the GMC 
based HyperSizer analysis are shown below. 



If this plot is compared to the non-GMC HyperSizer results or the Multiplies results 
shown in the previous section, the results are nearly identical. 

The strain results, which also match up very well with previous results are shown on the 
following page. 
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Again, this sample problem indicates that the HyperSizer implementation seems to work 
and appears to be feasible for even more large scale and complex panel concepts. The 
same argmuments hold for this case as far as checking the solution for accuracy by 
checking the upper and lower bounds on the solution using the total stiffness and elastic 
stiffness respectively. 

While this procedure has shown that the GMC implementation is feasible for a full panel, 
it also points out that there is still a considerable amount of work to be done to make it 
into a commerciable version. First, there is the disconnect between the stiffnesses used 
by the GMC and those specified in the *.VMT file. (For this case, the difficulty was 
overcome by making each ply homogeneous and making sure that the properties matched 
between the GMC file and the VMT file.) Second, the amount of data that must be 
generated for each ply in the GMC file is prohibitive and must be generated by hand. For 
this simple case of a single panel with only three plies, this is no problem, but for a 
problem with hundreds of panels, each with dozens of plies, the amount of input data 
would become unmanageable. 
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New Technology 


Several new and innovative techniques were developed and improved during the course 
of this research. 


First, the GMC method, originally developed by Dr. Jacob Aboudi, was extended by Dr. 
Carlos Orozco of UVA to include viscoelastic material analysis. In addition to just 
adding the capability. Dr. Orozco was able to substantially speed up the process and 
reduce the very large initial memory requirements. This was in an effort to make the 
code feasible for inclusion into HyperSizer. This capability should prove to be very 
valuable in performing very accurate, detailed analysis for fiber matrix composites. 
While we believe that there is still work to be done in order to make this capability 
commerciable, there were substantial advances and improvements made in this 
technology. 


Second, a substantial amount of effort went into the addition of the time-stepping 
algorithm into HyperSizer. While time stepping and iterative schemes are nothing new, 
their inclusion into HyperSizer is. This advance in HyperSizer went a long way into 
showing that this kind of analysis is completely compatible with HyperSizer's unique 
methods for analyzing and optimizing complex built-up stiffened panel and beam 
structures. In addition to adding the time-stepping to HyperSizer, a way of increasing the 
accuracy of the analysis without sacrificing accuracy was devised in the unique sub- 
iteration technique discussed in Appendix A. 

Finally, this research has demonstrated that viscoelasticity and GMC are feasible 
concepts for inclusion into HyperSizer. Furthermore, the methods developed open the 
door for more advanced technologies such as viscoplasticity and post-buckling strength 
analysis. While more work needs to be done to make these things commerciable, this 
research should prove to be a good base from which to work. 
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INCORPORATION OF A REVERSIBLE HEREDITARY VISCOELASTIC 

MODEL WITHIN GMC 

Carlos E. Orozco 

Department of Civil Engineering & Applied Mechanics 
University of Virginia, Charlottesville, VA 22903 


The following is a summary of the main relationships that I have used to incorporate the 
viscoelastic constitutive model of Ref. [1] within the context of the generalized method 
of cells (GMC). 


E m 

r~ 

Figure 1: Standard linear three element model. 


Equations (4) and (5) of Ref. [1] axe reproduced here in matrix form as: 



i a = c; 1 & a + e a t 

(1) 


e m = C;V m + 0 m T + £ v 

(2) 

where C 4 is the elastic stiffness matrix of a material with modulus of elasticity E a \ 
is the elastic stiffness matrix of a material with modulus of elasticity E m \ and the 8s 

c m 

are 

as follows: 

„ dc; 1 . 

«■ = gT o-. + w 

(3) 


acr 1 

dm - qj, O’m + w 

(4) 

where 

W = C‘ 3 “, “^," t 3 a \0,0,0} 1 ' 

(5) 

with 

™tan = w + Qf^ T 

(6) 
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In equations (1) through (6) subindex s refers to the single spring element in Figure 
1 and subindex m refers to the Maxwell element in Figure 1. w represents the coefficient 
of thermal expansion of the material. 

Solving for the stress rates in (1) and (2) we get: 

&, = C a (e, - 9 a T) (7) 

<rm = C m (e m - c v - 6 m T) (8) 

Now, since 


d* — <T m "f (f, 

Eqns. (7) and (8) can be combined to yield: 


(9) 


& = C E e - C m e v - C m 0 m T - C a 9 a T (10) 

where C E = C, + C m . 

Now, using the equal poisson ratio assumption of Ref. [1], the stiffness matrices can 
be written as: 

C m = E m N 

C a = E a N (11) 

C e = EN 

where E = E a + E m . 

Using (3) and (4), Eqn. (10) can be rewritten as: 


& = C E c — C m e v — C a 0T 


where 


( 12 ) 



0 = 0 + 9 — 1 rr 1 ^“'m 1 . 

s+ E a 6m dT 3 F E a dT m + 

E 

T, w 

(13) 

(Note that this is not the same 0 of Ref. [1].) 
Eqn. (12) can also be written: 




II 

0 

"fV 

1 

<v> 

-a 

1 

<V 


(14) 

where 

? Em . 

£ n- £ £v 


(15) 

and 

rp Wj - A • 

e T = =±0T 

E 


(16) 


Eqn. (14) constitutes a constitutive law that is analogous to that of plasticity or 
viscoplasticity. This form lends itself to easy implementation within the context of GMC 
without having to change the assembly of the matrices associated with the GMC proce- 
dure. 
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1 Formulation Without Internal Variables 

From equilibrium in the Maxwell element in Figure 1, we have: 


^ m — 7 l e v 

where, in analogy with (11) 

77 = 77 N 

Eqn. (12) can now be written as: 

& = Ce c - C m »7 “V m - C a 9T 

or, in view of (9) . . 

& — C E e - C m rj~ l ((r - tr a ) - C S 6T 

Solving now for e in (21), we get: 

. _ E m . . Em i ( x , 


c = %c -V + - *.) + 


which can be shown to be equivalent to Eqn. (17) of Ref. [1]. 


2 Solution of the Differential Equation 


Making use of the fact that: 


<r a = C.e 


and using relationships (11) and (19), Eqn. (21) can be written: 

& = C E i + — C,c C M (24) 

V V 

This is a first order differential equation that can be solved either for er or for e. 

A simple forward Euler scheme can be used to solve (24) by writing: 

<r(f + At) » «■(<) + |c E i(() + ^-C.e(() - ^V(i) - C,«(<)T j At (25) 
or, in view of (16), 

<r(( + A<) « o-(i) + |c E (i(i) - i T (t)] + ^C.s(t) - ^<r(*)} A t (26) 

where is the thermal strain rate. 

A sligthly more accurate way of solving Eqn. (24) numerically consists of finding 
its analytical solution assuming that either the temperature and strain history, or the 
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temperature and stress history is known. In what follows this is done for the former case 
(i.e., the temperature and strain histories are known). 

To this end rewrite Eqn. (24) as: 


where 


This is an O.D.E. of the form 
whose exact solution is: 


& + 7<r = ae + /?e — 6(t)T 

= En 

V 

a = Ce 

»=—c, 

V 

6(t) = C ,0{t) 
y' + f{x)y = r(x) 
y{x) = e~ h [J e h r{x)dx + C 


(27) 


where h = f f(x)dx. 

Therefore the solution of (26) is: 


<r(t) = e 7t jf ae^cdr + jT /?e 7T cdr — J 6(r)e' l ' T Tdr 


where 


m 


— c / gsi 

_c, \ dr 


E^d C- 1 £ 

J + E, dT m + E: 


: w } 


(28) 

(29) 


After some algebraic manipulations, an incremental version of Eqn. (28) is found as: 


<t(( + A() = e 7< “<r(<) + e (oi(i + ~) + + ^) - <(< + y)T(( + y)} Ai 

^ . < 3 °) 
This expression must be used instead of Eqn. (28) when the strain and temperature 

histories are not available in analytical form. 

In terms of the original parameters of the problem, Eqn. (30) can be rewritten as: 


<r(t + At) = e + e 


, A± 
2 


C E (c - i T ) + — C s c 
V 


At 


t+At/2 


(31) 


Note the similarity between Eqn. (31) and Eqn. (40) in Ref. (1]. This expression has 
been programmed whithin the new version of GMC for both the isothermal and the 
nonisothermal case. 

The figure that follows illustrates the results obtained with the new viscoelastic version 
of GMC. It corresponds to a relaxation test for TIMETAL 21S at 565 degrees C. It can 
be seen that the GMC approximation is very close to the exact solution. 

The nonisothermal case requires an additional approximation since the quantity S in 
(30) is not known at time t + At. This quantity is then approximated by <$(<). 
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3 Nonisothermal Case 


The nonisothermal case requires additional work on Eqn. (29) to make it suitable for the 
Euler solution procedure. 

Using (16) and (13) e T in (31) can be written as: 


; T = (§«. + ^-e m )T 


(32) 


where 9 a and 9 m are given by (3) and (4). 

Now, for a homogeneous isotropic material, Cj 1 corresponds to a compliance matrix 
given by: 

c ; 1 = Ej 1 L 

where 


L = 


1 —v —v 0 

—v 1 —v 0 

—v —v 1 0 

0 0 0 2 + 2i/ 


0 

0 


0 

0 


0 

0 


0 

0 


0 

0 

0 

0 

2 + 2i/ 

0 


0 

0 

0 

0 

0 

2 + 2i/ 


The partial derivative in Eqn. (3) can then be expressed as: 


where 

and 


dc; 1 _ j, F -x9L 

8T “ Ea dr L + La dT 


dL_du 

df~df Lo 


0 - 1-1000 
-1 0 -1000 
-1 -1 0 0 0 0 

0 0 0 2 0 0 

0 0 0 0 2 0 

0 0 0 0 0 2 


L 0 


The parameter 9 , can now be expressed as: 


9 a = (- E ; 2 ~r L + E; x ^L LoK + w 


The analogous expression for 9 m is then: 

,dE 


Om = (-£, 


'm 


-2 

m Qj 


du 

L + £-‘gj;L„)<T m + W 


(33) 


(34) 


(35) 

(36) 


(37) 


(38) 


(39) 


It is necessary now to find expressions for <r 3 and <r m to complete the derivation. 
Given the stress rate <T and the thermal strain rate e T at a given time t, the new modified 
viscous strain can be obtained from (14) as: 


«r, = -C E l & + e — c 
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(40) 



Once the new total stress and viscous strain rate are known, tr 3 and tr m 
obtained as: 


E „ t 

r m — ^ CtjC 


E„ 


V'-'V 


can then be 
(41) 


and 


= -a m 


(42) 


The preceding derivation provides a complete account of all the expressions needed to 
implement an Euler-based numerical integration procedure for the multiaxial nonisother- 
mal viscoelastic model as presented in [1], 
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Progress on laminate/panel viscoelastic integration 
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Abstract 

Task number 2 of the NRA contract involves incorporating fiber/matrix composite 
viscoelastic behavior into the commercial sizing/analysis code called HyperSizer. 
Dr. Carlos Orozco has advanced this research by applying the General Method of 
Cells (GMC), originally developed by Dr. Jacob Aboudi, to the viscoelastic 
behavior of general composite materials. In the current effort, an algorithm to 
incorporate this GMC methodology with either a single ply, a multi-ply laminate 
or a full stiffened panel has been developed and is currently under development 
and testing. Because the analysis of a full structural panel is a complex problem 
with a limited number of available exact solutions, a simplified formulation was 
developed to predict and analyze the viscoelastic behavior of 1-D strip elements 
composed of multiple “plies.” By analyzing a simplified problem for which 
closed form, analytical solutions were known, the numerical algorithm was tested 
and its performance was characterized, and is presented herein. The numerical 
algorithm for a full panel subjected to specified constant stress (“creep”), 
specified constant strain (“stress relaxation”), specified transient stress or specified 
transient strain loading is presented. A number of closed form solutions are 
presented against which the simplified numerical formulation was tested in the 
form of a stand-alone, prototype computer code. The comparisons indicate that 
the methodology is consistent and 1 st order accurate, that is, the solution accuracy 
is directly proportional to the magnitude of the time step. Also, while the forward 
Euler time stepping scheme is conditionally stable, the time step restriction is 
easily characterized as a function of the viscoelastic material properties and can be 
easily automated. In addition, a sub-iteration scheme was developed to allow 
multiple time steps to be performed for the individual plies before advancing the 
global time step. This minimizes the number of switches between different 
materials and plies. This modification to the formulation is shown to improve the 
efficiency of the analysis without compromising the accuracy. The next steps are 
to include the GMC subroutine and full panel formulation into the prototype code. 
This will allow characterization of GMC behavior with the numerical formulation 
presented here. The final step is to integrate GMC and the viscoelastic 
methodology into HyperSizer itself. 
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1. Introduction 

One of the primary tasks of the NRA contract was to incorporate the effects of 
viscoelasticity into the commercial analysis / optimization code called HyperSizer. 
In particular, we wish to model the viscoelastic behavior of orthotropic materials 
such as polymeric or metal matrix composites. In order to accomplish this, we are 
using the General Method of Cells (GMC) originally developed by Dr. Jacob 
Aboudi to model the micromechanical behavior of the fiber/matrix composites. 
This research is aimed at incorporating the micromechanical response obtained 
using GMC into a macromechanical model predicting the behavior of a single ply, 
a laminate consisting of multiple plies, a panel consisting of multiple laminates 
(such as a corrugated panel) and ultimately, a structure composed of multiple 
panels. 

The first step in this integration lies in determining the best way to use the GMC 
method to extend from a single ply to a laminate composed of multiple plies. In 
order to accomplish this, it is useful to take a step back from the full plate theory 
analysis to a simpler 1-D analysis involving a simple bar exhibiting “strip” 
behavior as opposed to “plate” behavior. This bar could be composed of one or 
many plies, exposed to either tension or compression. This simple analysis allows 
us to derive analytical or “exact” solutions to which we can compare our 
numerical methods for analyzing laminates and thus determining if the method 
that we are using is correct. 

This report details the development of an algorithm for integrating the 
micromechanical, single material analysis methodology into a laminate response. 

A series of analytical solutions were derived and compared to the simplified 1-D- 
numerical model of this algorithm in a computer code and the results are presented 
and discussed below. The final section contains a discussion of future work and 
the next direction for this research effort. 

Before going into detail about the work performed here, we present a short 
background leading to the development of the proposed algorithms. 

2. Background 

Viscoelasticity is a complex time-dependent phenomena that becomes very 
difficult to model in closed form for any type of composite material. Dr. Carlos 
Orozco has extended the GMC methodology of Dr. Aboudi to perform viscoelastic 
analyses using a simple forward Euler scheme to step the viscoelastic solution in 
time. His implementation is based on an observation that the stress rate, a in a 
viscoelastic member can be expressed as: 
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& = e{s + S n + E T ) 



where E is the overall stiffness, E m is the viscoelastic stiffness, e is the overall 
strain rate, s n is the strain rate in the viscous damper, and s T is the thermal strain 

rate. The period superscript notation (i.e. e ) indicates derivative with respect to 
time ( ds/dt ). When cast in this form, it is apparent that the viscous strain 

quantity, acts in exactly the same manner as the thermal strain, therefore, 

Carlos’ idea is to implement the viscous strain for a laminate in the same well 
established way that HyperSizer implements thermal strain. 

Carlos’ implementation of GMC takes the form of a FORTRAN subroutine which, 
at each time step, returns the quantity = f{At,e n ,e„^) where the subscript n 

denotes the time step number and At is the magnitude of the current time step. 

The quantity e n , calculated for each ply, is then assembled into a laminate 

viscoelastic response using an equation similar to that given in equation 10 from 
Ref. [1] This equation is used to assemble ply level thermal strains into a laminate 
level thermal response. The proposed algorithm is listed in detail in the 
Formulation section below. 

3. Purpose 

The purpose of the current effort is to: 

• Develop an algorithm for integrating the GMC micromechanical viscoelastic 
model into the HyperSizer commercial sizing/analysis software. 

• Demonstrate the effectiveness of this algorithm by developing a simplified 1-D 
analog and applying the simplified model to problems with known solutions 
involving 1 or 2 viscoelastic or elastic plies. 

• Develop a testbed that can be used to test modifications to the proposed 
algorithm which may improve the efficiency. For example, one such 
modification is the addition of an arbitrary number of sub-steps to be 
performed at the ply level before updating the laminate. With the current 
implementation of the GMC subroutine, this has the advantage of limiting the 
number of times the GMC code must switch from one ply material to another 
and thus limiting the number of matrix inversions. The simplified prototype 
model gives a way of testing such proposed modifications. 

• Provide a framework for implementing an isotropic viscoelastic material in 
HyperSizer independent of the GMC procedure. This will involve extending 
the simplified 1-D method to include aspects of plate and/or beam theory such 
as bending and Poisson effects. 
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4. Formulation 

There are two algorithm formulations presented. For each formulation, there are 
two sub-cases. The first algorithm is for an implementation of a full panel 
formulation as in HyperSizer. The second formulation involves the simplified 1 -D 
analog from which the exact results are taken. For each formulation, the sub-cases 
involve either specified stress fields or specified strain fields. 

In addition to the algorithm formulation, there were a number of analytical or 
‘exact” solutions identified for comparison to the numerical algorithms. These 
exact solutions are identified and listed in appendix A. 

4. 1 HyperSizer (full panel) Formulation 

The full panel formulation is divided into the sub-cases of either response to a 
specified stress field or specified strain field. 

4.1.1 Specified Stress Formulation 

In this scenario, the stress field (or panel forces) is specified as a function of time 
and the algorithm determines the resulting strain field. A special case of this type 
of analysism called “creep,” occurs when the specified stress field is constant. The 
algorithm is laid out in the attached Flowchart 1. 

In this (and all following flowcharts), the subscript capital N represents the global 
time step, lower case n represents the time step for each sub-iteration. All 
quantities with lower case k subscripts represent single ply level quantities and 
quantities without this subscript are laminate level quantities, (i.e. } is a single 

ply strain field while {£•} represents a laminate strain field) 

The following notes apply to the step numbers shown in Flowchart 1 : 

1 • The first step is to calculate the elastic response of the laminate given the 
change in external forces from the previous to the current time step. [C] is the 
laminate stiffness matrix. 

2. Using the strain field at the current time step, calculate the individual ply 
strains using the laminate strain and curvature field. 
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Flowchart 1 : HyperSizer (full panel) specified stress viscoelastic formulation 
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3. Using either GMC or some form of standalone subroutine, calculate the 
viscous strain quantity e n . This step is performed at each sub-iteration and the 

time step passed to this function is the time step associated with the sub- 
iteration (i.e. At n = At / (# of sub-iterations)). The rationale for using these sub- 
iterations is to allow multiple time steps for each ply while limiting the number 
of times a new ply must be loaded. For a homogeneous viscoelastic material, 
the GMC routine is not required to advance the time step, and the viscous 
strain quantity is calculated from: 

trj 

4. Calculate the change in the viscous force from the previous global time step to 
the current global time step. 

5. This equation is similar to Equation 10 from Reference [1]. This equation 
allows us to “assemble” the viscous strain contributions from each ply into a 
laminate viscous force. This is exactly how ply level thermal strains are 
currently “assembled” into laminate thermal forces and moments. 

6. Update the laminate strains to reflect the change in the viscous forces. Note 
that the viscous strain calculated in step 3 is dependent on both the current and 
the previous levels of strain. Therefore, it may be desirable to feed these new 
updated strains back to step 2 and iterate 1 or more times. It is unclear at this 
time whether this iteration will be necessary. 

7. Advance to the next time step and return to step 1 . 

4.1.2 Specified Strain Formulation 

In this scenario, the strain field is specified as a function of time and the algorithm 
determines the resulting stress (or force) field. A special case of this type of 
analysis occurs when the specified strain field is constant. This situation is called 
“stress relaxation”. The algorithm is laid out in Flowchart 2. 

The algorithm is very similar to the algorithm for specified stress with the 
exception that instead of updating the strain field at each time step, this algorithm 
aims to update each plies stress and resulting laminate internal forces and 
moments. Note that the loop performed at the ply level (steps 3 through 5) which 
calculates the viscous forces is exactly the same as that of the specified stress 
formulation. 
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Flowchart 2: HyperSizer (full panel) specified strain viscoelastic formulation 
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Flowchart 3 : Simplified model specified stress (creep) viscoelastic formulation 
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4.2 Simplified Formulation 

The simplified formulation attempts to model a problem involving 1-D 
deformation of a laminate composed of multiple plies in the same way as the full 
panel formulation. In this case, 3-D effects including curvature and Poisson’s 
effect are neglected and the material is characterized by a single physical stiffness, 
E. If the ply is viscoelastic, the stiffness, E, which we can think of as the 
“instantaneous” stiffness is the sum of the elastic stiffness, E s plus the viscoelastic 
stiffness, E m . Also note that for this analysis, there is only one strain as each ply 
strains at the same rate as the laminate, therefore it makes no sense in the 
simplified analysis to speak of laminate vs. ply strains. Note that in the simplified 
analysis, the laminate “stiffness matrix” reduces to the physical stiffness quantity, 
E, multiplied by the laminate thickness, H. In addition, the integral in step 5 (Eq. 
10 from Ref. [1]) reduces to: 


U plies 


(AW')= 


k = 1 


4.2.1 Specified Stress Formulation 

A special case for this formulation occurs when the laminate stress is constant. 
This special case is called “creep.” The formulation, laid out in flowchart number 
3, is very similar to the specified stress formulation given for the full panel above. 
One special note about the simplified algorithm. By using the forward Euler 
method, the procedure in step 3 to calculate the viscous strain, e , is actually 

independent of the strain at the current time step. In other words, the outer loop 
(after step 6) is unnecessary for the simplified algorithm, however the step is left 
in the flowchart to maintain similarity with the full panel formulation. 

4.2.2 Specified Strain Formulation 

This formulation, laid out in flowchart number 4, is exactly analogous to the 
specified strain formulation for the full panel above. 

4.3 Analytical Results for Simplified Analysis 

The following matrix illustrates the cases for which analytical or “exact” results 
were derived for comparison to the numerical results from the previous section. 
There are nine basic cases involving either one or two plies. 
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Table 1 : Loading cases for comparison of analytical and numerical results 



Creep 

Analysis 

Stress Relaxation 
Analysis 

Transient Strain 
Profile 

1 Viscoelastic Ply 

Analytical 

Analytical 

Analytical 

1 Viscoelastic & 
1 Elastic Ply 

Analytical 

Analytical 

Analytical 

2 Viscoelastic plies 

Runga-Kutta 4th Order 

Analytical 

Analytical 


The analytical solutions are listed in detail in Appendix A. Note that for the case 
of creep analysis in 2 viscoelastic plies, instead of solving the resulting system of 
2 differential equations in closed form, we chose to model this response using a 
4th order Runga-Kutta analysis which, while not exact, is much more accurate 
than the 1st order Euler numerical integration used in the present effort. 

5. Results / Discussion 

The simplified formulation was implemented in a Windows computer program 
that simulates any number of plies and any combination of viscoelastic or elastic 
plies with arbitrary thicknesses. The ply materials used in the analysis are listed in 
Figure 1 . The first material is a generic fiberglass, modeled here as an elastic 
material with a stiffness of 68, 940 MPa. The program is written to be general but 
no provision is made to model a purely elastic material. Therefore, in order to 
model an “elastic” material, the viscoelastic damping parameter was set to a very 
large value (1 .Ox 10 20 ). While this is technically a viscoelastic material, the 
analysis would have to be run for approximately 10 15 seconds in order to notice • 
any viscoelastic response. Therefore, for all practical purposes, the material is 
elastic. The second and third plies are both composed of a viscoelastic epoxy 
taken from Dr. Aboudi’s book. The difference being in the material evaluation 

temperatures of 22°C (cold) and 140°C (hot). 

Figures 2-1 1 are result plots using combinations of these three plies under the 
loading conditions laid out in Table 1. In each of the plots, the circle makers 
represent the numerical results while the solid lines represent the exact or 
analytical result for comparison. The “maximum error” shown for each result is 
calculated from: 


C nux = MAX 

P ~P 

x numerical exact 

P 


exact 


The obtained results demonstrate that the method is consistent, that is, it converges 
on the exact solution as the time step is reduced to zero. Figures 2a-e illustrate 
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this trend as well as the effectiveness of the sub-iteration technique. Because the 
chosen method of time integration (forward Euler) is first order accurate, one 
would expect the error to be linearly proportional to the size of the time step. In 
figures 2a and 2b, the ratio of the time steps is 5.0, and the ratio of the maximum 
errors is about 5.8. Reducing the time step by another factor of 5.0 in figure 2e 
reduces the maximum error by a factor of 5.2. This behavior continues as the time 
step size is reduced to zero, as illustrated in Figure 2g. 

a) Elastic ply 


Epoxy (cold) New Ply 

|3 4470 Elastic StiKness 

tpoxy (Mot) 

DeMaPty 

[34470 VtscoSWfn«S8 


1 OOOE+20 Damping Coefficient 

jo. 01 000 Thickness 

i |2 901E*15 Tnw COMtant 

Close 

[68940. SMuw 


b) Cold viscoelastic ply 


Epo>y(Hot) 


New Ply 
Delete Ply 


Close 


c) Hot viscoelastic ply 


Gloss-V 
El 


New Ply 


Delete PV 


Qose 


J4057. 

Elastic Sfcftwt* 

[199.0 

VIscoS teteest 

|l 47000 

Damping CoafStiant 

[O.OSOOO ~ 

Thickness 

j?38.7 

““ Time Constant : ; 

|4256. 

ovinviv 


(3080 ~~ 

rmtr«im» ~ 

(335.0 

“ ViscoStitees* 

(400000. 

Damping Coefficient 

[0.05000 

Thickness 

(l194. 

Time Constant 

(3415. 

" Stifoees 


Figure 1: Ply properties using in viscoelastic analysis tests 
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In addition to the convergence of the algorithm with 
decreasing time step, the sub-iterations are also effective 
in reducing the error. As shown in Figures 2a and 2b, 
reducing the time step by 5 decreases the error by a factor 
of 5, but Figure 2c demonstrates that increasing the 
number of sub-steps by a factor of 5 has almost the same 
effect. In fact the reduction of error obtained in this case 
is more than 1 0 times, meaning that this is a better 
solution than that of Figure 2b. While the results 
obtained by increasing the number of sub-iterations are 
not always better, it was observed that they were always 
comparable. It appears that there is some “optimum” 
value for the number of sub-iterations which we did not 

determine, but between 5 and 10 generally gave very good results. One word of 
caution. One of the drawbacks of using Euler integration is that the method is 
conditionally stable. If the time step is too large, oscillations in the solution occur 
as illustrated in Figure 2f. The practical limit on the time step size is the minimum 
time constant of the viscoelastic materials in the laminate. This time constant, x, is 
calculated from the viscoelastic properties as: 


All Figures are 
sealed to show 
the normalized 
percentage ol 
change to quickly 
assess the relative 
significance ol 
the viscoelastic 
behavior. 


T = 


JL 

E_ 


where rj is the viscoelastic damping coefficient and E m is the viscoelastic stiffness. 
For the material shown in Figure 2f, the time constant is 738 seconds and the time 
step used is 1000 seconds. 


Although convergence results are not shown for the rest of the cases, all of the 
loadings and ply combinations exhibit the same convergence and stability 
behavior as that shown in figure 2. 


Figures 2-1 1 demonstrate 
the solution accuracy. 

The numerical solution is 
identified with markers 
while the exact solution is 
identified w ith a solid 
line. 


Figure 3 shows the stress relaxation result for a 
single viscoelastic ply. Figures 4 and 5 illustrate 
a single viscoelastic ply subjected to a specified 
transient strain (the strains are plotted in Figures 
4b and 5b respectively). Figure 4 shows the 
result as the ply is subjected to a 10% tensile 
periodic loading and Figure 5 illustrates the 
response under a 150% loading where the 
loading oscillates between tension and 
compression. 
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Figure 2a: Creep Analysis with 1 viscoelastic ply (At = 500 s; 1 sub-iteration) 



Figure 2b: Creep Analysis with 1 viscoelastic ply (At = 100 s; 1 sub-iteration) 
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Figure 2f: Oscillations resulting from large timestep 


NAS A/CR— 2001-21 1 166 


38 









log(Max. Error) 



Figure 2g: Convergence of Creep result for a single viscoelastic ply 



Figure 3: Stress relaxation result with 1 viscoelastic ply (At - 100 s; 5 sub- iterations) 
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Time step 











Figure 5a: Specified transient strain solution (larger magnitude At = 100 s; 5 sub-iterations) 



Figure 5b: Transient strain profile for analysis 
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Figures 6 and 7 illustrate creep results for the case of a viscoelastic and elastic ply 
and the case of 2 viscoelastic plies respectively. Because no closed form solution 
was obtained for the case of creep in two viscoelastic plies, the result was 
compared to a result obtained using a 4 th order Runga-Kutta integration technique. 
As the 4 th order technique is assumed to be much more accurate them Euler (which 
is 1 st order accurate), the Runga-Kutta solution was assumed to be sufficient for 
comparison. 

Figures 8 and 9 illustrate stress relaxation response for the case of a viscoelastic 
ply and elastic ply and the case of 2 viscoelastic plies respectively. Figures 10 and 
1 1 demonstrate the same two cases under a specified transient strain loading 
profile. 

Figures 12-15 illustrate some interesting features of the viscoelastic analyses 
involving multiple plies. In Figures 12 and 13 a laminate composed of one elastic 
and two viscoelastic plies is subjected to creep (constant stress) and stress 
relaxation (constant strain) loading respectively. In the case of creep, one can see 
that over time, the stress in the viscoelastic plies decreases as more and more of 
the stress is picked up by the elastic ply. Also, notice that the hot ply, which 
exhibits a stronger viscoelastic behavior picks up less of the total load than that of 
the cold viscoelastic ply. In Figure 13, the stress in the hot ply relaxes more than 
that of the cold ply, and of course, the stress in the elastic ply is constant. Figures 
14 and 15 illustrate the same load cases with a laminate made up of two 
viscoelastic plies. Figure 14 is interesting because it shows that over time, the 
cold epoxy ply (which is less viscoelastic) exhibits an increase in stress while the 
hot epoxy experiences a decrease in stress. The net effect is a zero change in the 
laminate stress over time. 
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Figure 6: Creep Solution for one viscoelastic and one elastic ply (At = 100 s; 5 sub-iterations) 
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Figure 8: Stress relaxation result for one viscoelastic and one elastic ply (At = 100 s; 5 sub-iterations) 



Figure 9: Stress relaxation result for two viscoelastic plies (At = 100 s; 5 sub-iterations) 
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Figure 12: Individual ply stresses for two viscoelastic (epoxy) plies and one 
elastic ply (glass) undergoing creep 



Figure 13: Individual ply stresses for two viscoelastic (epoxy) plies and one 
elastic ply (glass) undergoing stress relaxation 
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6. Conclusions 

Several conclusions can be drawn from the preceding results: 

1. For the simplified problem of 1-D deformation in a laminate, the 
method of combining the responses of multiple plies to obtain a 
laminate response appears to be consistent. In other words, it converges 
on the correct result as the time step is reduced to zero. 

2. The method is conditionally stable. The restriction on the time step is 
approximately equal to the minimum viscoelastic time constant, x, of 
the materials in the laminate. 

3. The sub- iteration scheme, which allows several small time-steps on 
each individual ply before assembling those results into the laminate is 
effective in improving the overall accuracy while providing increased 
computational efficiency. 

7. Future Work 

The current work indicates that it is worthwhile to pursue the full panel 
formulation, eventually resulting in integration of viscoelastic plies into 
HyperSizer. There are two enhancements that are planned to be included in a 
separate "stand-alone” prototype code before attempting to integrate with 
HyperSizer. The first enhancement involves extending the present work to a full 
3-D laminate including the effects of laminate curvature and Poisson effects. This 
algorithm will look very similar to that for a full panel as given in flowcharts 1 and 
2. Second, it would be useful to integrate the Fortran subroutine, GMCVISC into 
the standalone code to determine the characteristics of using GMC in the 
algorithm before attempting to integrate with HyperSizer. 

It also appears that the viscoelastic issues addressed and solved under the NRA 
funding will be directly applicable to developing a future viscoplastic laminate and 
panel capability. 

The final phase of this effort is to integrate this algorithm into HyperSizer itself. 
The current implementation of the algorithm (with 3-D effects) can most likely be 
integrated directly into the code without GMC to model homogeneous viscoelastic 
materials (i.e. adhesive bond-coats, thermal barrier coatings, etc.) The inclusion of 
the viscoelastic algorithm with or without GMC will involve the inclusion of a 
time parameter and may be a substantial effort. At the present time, HyperSizer 
performs point analyses which are independent of time. A viscoelastic capability 
for polymer composite ‘built-up’ panels will be a substantial and valuable added 
capability in HyperSizer. 
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Appendix A Analytical Results for Viscoelastic Formulation 

A.1 Results for a Single Viscoelastic Ply 


A. 1.1 Governing Equation 

We are modeling a viscoelastic material using a 
modified Maxwell model as given in the Figure on the 
right. The quantity E s represents the viscoelastic 
stiffness, E m represents the viscoelastic stiffness and 7 
represents the viscoelastic damping. The basic 
equation governing the overall response of this 
material is: 



where £ is the strain rate and <y m is the stress in the viscoelastic stiffness. 

A. 1.2 Stress Relaxation 

For a constant strain value given by s 0 , the stress is given by: 

a{t)=£ Q (E s +E m e- r ') 

where y is the inverse of the material viscoelastic time constant and is given 
by: 


Y 



A. 1.3 Creep 

For a constant stress, cr 0 , the strain is given by: 


K')= — 

E. 


1 - 


E.+E. 


where in this case, the time parameter is: 


\ 

J 


y. E - E - 

(E. + E,)n 
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A. 1.4 Specified Transient Strain 

Assuming a given linear strain profile given by: 

e{t) =mt + e 0 

The stress as a function of time can be expressed as: 

cr(t) = e(t)E s + \e 0 E m e- Y ' + mr}{ 1 - )J 

where y=^E m /rj. 


A.2 Solutions for One Viscoelastic Ply and One Elastic Ply 


A.2.1 Governing Equations 

In this case, the system is modeled by the 
system shown at the right. The governing 
equations in each branch of the system are 
given by: 





A.2.2 Stress Relaxation 

Because the strain field is specified and always known, the viscoelastic stresses 
act independently of the elastic stresses, so the integrated stress is actually a 
linear combination of the stress in each branch of the model. Also, it makes 
more sense to speak of overall forces for the system since the two branches 
(elastic and viscoelastic) can have different areas (or thicknesses). The 
solution for force in this two ply laminate is given by: 

P = Pe+P, 

where y is the inverse of the time constant for the viscoelastic ply, y = E m /r). 
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A.2.3 Creep 

The solution for specified constant force, P 0 is given by: 


K')= 




A.E. + A.E, + A,E„(\-e- r 'f 
a,e, + a(e, + e.) _ 


where the time parameter is given by: 

e,(a.e,+a,e,) 


r = 


n (A,E,+A,(E,+E.)) 


A.2.4 Specified Transient Strain 

As in the case of stress relaxation, because the strain is known, the elastic and 
viscoelastic sides act independently and the solution for force is given by. 


P{t)=P e + P v 

P{t)= A e E e £(t)+A v [E s e{t)+£ 0 E m e- r ‘ 

where the time constant is taken from the viscoelastic branch as y = E m / ij . 


A.3 Results for Two Viscoelastic Plies 


A. 3.1 Governing Equations 

In this case, the system is modeled by the 
system shown at the right. The governing 
equations in each branch of the system are 
given by: 


£ 


Vl 







A.3.2 Stress Relaxation 

As in the stress relaxation case outlined in A.2.2, the integrated stress for this 
two ply laminate is a linear combination of the independent stresses in each 
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branch of the model. The solution, cast in terms of the laminate internal force, 
is given by: 


p(i)=p,+p ! 

/>(<)= £„ U k, + )+ a, (E n + )] 

where y, = and y, = P.Jn, 


A.3.3 Creep 

The case of creep in a laminate composed of two viscoelastic plies is the most 
complex solution presented here. The governing equations applied to the 
problem of creep reduce to: 


£■(/) = — fjnl !jo2 
A\E S \ + A 2 E s2 

where the P 0 is the complete force applied to the laminate and the forces in the 
viscoelastic branches ( P ml and P m2 ) of each ply are determined from the system 
of first order differential equations: 


where: 


dP ml ] 


-c, 

d \C 2 

dt 


0-AA) 

(i-AA) 

dP m2 


A c, 

-A 

l dt \ 


„0-AA) 

(l-AA) 


q _ E m \ + E^Aj) 

Q _ + &s2^2 ) 

r h{Es 2 A 2 +A l {E s] +E ml )) 


D. - A ' E ” X 

jT) _ ^2^ m2 

E s 2 A 2 + Aj (E sl +E m] ) 

‘ £„4 + A, (£„+£.,) 


This set of differential equations was not solved in closed form in the present 
work. Instead, a 4 th order Runga-Kutta numerical integration scheme was used 
to find a solution which, while inexact, is much more accurate than the 
numerical solution resulting from the 1 st order forward Euler algorithms 
presented here. The actual Runga-Kutta implementation is beyond the scope 
of this report, but can be found in any text on numerical methods, (e.g. Ref.[2]) 
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A.3.4 Specified Transient Strain 

As in the case of stress relaxation, the two branches of the model act 
independently and the solution is given by: 

p if)= p ^ p i , V. 

P(t)= A l [E s] e(t)+e 0 E ml e- r '‘ A 2 [E s2 e{t)+e Q E ml e- ri ‘ + w^ 2 (l-e' rj ')J 

where the time constants are y x = E mX jrj x and y 2 = • 
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